#ifndef __DISCRETEFORM__
#define __DISCRETEFORM__

#include "femspace.h"
#include "femfunction.h"
#include "mesh.h"
#include "base.h"
#include "enumerations.h"
#include "constants.h"
#include "quadrature.h"

typedef struct AUXFEMFUNCTION_
{
    /** 需要多少个有限元函数 */
    INT NAuxFEMFun;
    /** 存储有限元函数列，可能会包含多个有限元函数 */
    FEMFUNCTION **AuxFEMFun;
    /** the number of Multiindex for the AuxFefunctions */
    /** 对这些有限元函数我们需要对他们取什么样的导数值*/
    /** 就是说AuxFEMFun中的每个有限元函数需要在数组AuxFEMFunMultiIndex中需要计算的导数个数 */
    // 比如 NAuxFEMFunMultiIndex=[3，2，1]: 表示的是前3个导数值是第0个有限元函数的导数
    //  接下来的2个导数值是第1个有限元函数的导数, 最后一个导数是第2个有限元函数的导数
    INT *NAuxFEMFunMultiIndex;
    /** 对有限元函数需要取值的导数指标 */
    // 记录的是一组导数的符号：比如[D00 D10 D01 D11 D00 D10 D01]
    MULTIINDEX *AuxFEMFunMultiIndex;
    /** 从有限元函数列中取到的在某一点的函数值的个数 */
    /**这个应该可以直接计算出来的，不需要给定：*/
    // 具体的值是：NAuxFEMFunValues = sum_i(NAuxFEMFunMultiIndex[i]*dim_values[i]);
    INT NAuxFEMFunValues;
    // 下面这个实数数组用来存储辅助有限元函数计算值的时候所需要的内存空间
    DOUBLE *AuxBaseValues;
    /** 用来存储有限元函数值的数组： 维数 NAuxFEMFunValues */
    DOUBLE *AuxFEMFunValues;
    // lhc 添加
    // 是否存在提前计算出的标准单元积分点上的信息
    BOOL IsHaveAuxBaseValuesInPointsOfRefElem;
    // 下面三个信息的指示子 表示MultiIndex中1，2，3阶个数和总个数
    INT IndAuxMultiIndex[4];
    /** 如果基函数只依赖于RefCoord 建立储存RefElem中积分点上有限元基函数的所有微分算子的值, 长度是 num_point*num_base*valuedim*NumMultiIndex */
    DOUBLE *AuxBaseValuesInPointsOfRefElem;
    // 积分点编号 由于ComputeAuxFEMFunValue函数中不包含积分点信息 无法确定从AuxBaseValuesInPointsOfRefElem中读取的位置
    INT Points_ind;
} AUXFEMFUNCTION;

AUXFEMFUNCTION *AuxFEMFunctionBuild(INT NFEMFuns, FEMFUNCTION **femfunctions,
                                    INT *NAuxFEMFunMultiIndex, MULTIINDEX *AuxFEMFunMultiIndex);
void AuxFEMFunctionDestroy(AUXFEMFUNCTION **AuxFEMFunction);
// 计算辅助函数在积分点的函数值
void ComputeAuxFEMFunValue(AUXFEMFUNCTION *AuxFEMFun, INT ElemInd, ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[]);
// 计算多个有限元函数之间的运算的值: 比如$\|u_{1,h}-u_{2,h}\|_{1,\Omega}$
// 我们建设所有的有限元空间是定义在同一个网格上
DOUBLE ComputeFEMFunctionValue(AUXFEMFUNCTION *AuxFEMFun, ERRORFORM *ErrorForm, QUADRATURE *Quadrature);

typedef struct DISCRETEFORM_
{
    /** the anasatz space, corresponding to the column index */
    /** 取值有限元函数空间，对应矩阵的列 */
    FEMSPACE *LeftSpace;
    /** the Multiindex for testspace */
    /** 形成刚度矩阵的时候需要检验有限元空间的导数个数 */
    INT NumLeftMultiIndex;
    /** 形成刚度矩阵的时候需要检验有限元空间的导数 */
    MULTIINDEX *LeftMultiIndex;
    /** 建立存储取值有限元基函数的所有微分算子的值, 长度是 num_base*valuedim*NumLeftMultiIndex */
    DOUBLE *LeftBaseValues;
    /** test space， corresponding to the row index */
    /** 检验有限元函数空间，对应矩阵的行 */
    FEMSPACE *RightSpace;
    /** the Multiindex for anasatzspace */
    /** 形成刚度矩阵的时候需要试探有限元空间的导数个数 */
    INT NumRightMultiIndex;
    /** 形成刚度矩阵的时候需要试探有限元空间的导数 */
    MULTIINDEX *RightMultiIndex;
    /** 建立存储试探有限元基函数的所有微分算子的值, 长度是 num_base*valuedim*NumLeftMultiIndex */
    DOUBLE *RightBaseValues;
    /**建立边界条件函数　*/
    FUNCTIONVEC *BoundaryFunction;
    // 230404 add
    FUNCTIONVEC *NeumannBoundaryFunction;
    FUNCTIONVEC *RobinBoundaryFunction;
    /** 处理Neumann和Robin边界条件时需要的边界上的积分格式 */
    QUADRATURE *PartialQuadrature;
    // 变分形式中需要的辅助有限元函数(服务于非线性问题)
    AUXFEMFUNCTION *AuxFEMFunction;

    // lhc 添加
    // 是否使用非线性对应的有限元网格进行组装
    BOOL IsUseAuxFEMSpace;
    // 是否存在提前计算出的标准单元积分点上的信息, Ciarlet三元组的时候可以提前计算标准单元上的导数值
    BOOL IsHaveLeftBaseValuesInPointsOfRefElem;
    BOOL IsHaveRightBaseValuesInPointsOfRefElem;
    // 利用非线性函数中的哪个有限元空间进行组装 如果用原来的左有限元空间则此为NULL
    FEMSPACE *LeftSpaceInAssemble;
    // 两个FEMSPACE的层级关系到底用哪个
    INT *Relates;
    // 下面三个信息的指示子 表示MultiIndex中1，2，3阶个数和总个数, 最后一个位置存储总共用了多少个基函数的微分值
    INT IndLeftMultiIndex[4];
    INT IndRightMultiIndex[4];
    /** 如果基函数只依赖于RefCoord 建立储存RefElem中积分点上有限元基函数的所有微分算子的值, 
     * 长度是 num_point*num_base*valuedim*NumMultiIndex */
    DOUBLE *LeftBaseValuesInPointsOfRefElem;
    DOUBLE *RightBaseValuesInPointsOfRefElem;

    /** Discrete Form for the Matrix */
    /** 定义形成刚度矩阵的变分形式 */
    DISCRETEFORMMATRIX *DiscreteFormMatrix;
    /** 定义行程右端项的变分形式 */
    DISCRETERHSFORMVEC *DiscreteFormVec;
    /** The information of the quadrature scheme */
    /** 形成刚度矩阵需要的积分点的信息 */
    QUADRATURE *Quadrature;
    // 下面的空间是专门用来组装刚度矩阵的时候用来做数据交换的内存空间
    DOUBLE *AuxValues;
    // 表示目前是否是计算特征值问题中的质量矩阵
    BOOL IsMassMatrix;
} DISCRETEFORM;

// 下面是一些操作函数：最重要的是建立离散变分形式
/** 建立离散变分形式 （不需有限元函数）*/
DISCRETEFORM *DiscreteFormBuild(FEMSPACE *LeftSpace, INT NumLeftMultiIndex, MULTIINDEX *LeftMultiIndex,
                                FEMSPACE *RightSpace, INT NumRightMultiIndex, MULTIINDEX *RightMultiIndex,
                                DISCRETEFORMMATRIX *DiscreteFormMatrix, DISCRETERHSFORMVEC *DiscreteFormVec,
                                FUNCTIONVEC *BoundaryFunction,
                                QUADRATURE *Quadrature);
/** 建立离散变分形式对象 （需要有限元函数）*/
DISCRETEFORM *DiscreteFormAuxFeFunBuild(FEMSPACE *LeftSpace, INT NumLeftMultiIndex,
                                        MULTIINDEX *LeftMultiIndex, FEMSPACE *RightSpace,
                                        INT NumRightMultiIndex, MULTIINDEX *RightMultiIndex,
                                        AUXFEMFUNCTION *AuxFEMFunction,
                                        DISCRETEFORMMATRIX *DiscreteFormMatrix, DISCRETERHSFORMVEC *DiscreteFormVec,
                                        FUNCTIONVEC *BoundaryFunction,
                                        QUADRATURE *Quadrature);
// 释放内存空间
void DiscreteFormDestroy(DISCRETEFORM **DiscreteForm);

// lhc添加
/* 提前计算一些单元公共的信息 省略重复计算 */
// 计算RefCoord中各个积分点对应的函数值信息(只计算能用于简化计算的部分)
void ComputeBaseValuesInPointsOfRefElem(DISCRETEFORM *DiscreteForm);
// 释放计算出的RefCoord中各个积分点对应的函数值信息
void DestroyBaseValuesInPointsOfRefElem(DISCRETEFORM *DiscreteForm);
/**
 * @brief 根据Elem中需要的导数信息获得RefElem中的导数信息
 * @param MultiIndexInElem 输入计算实际单元的函数值时需要的导数信息
 * @param NMultiIndexInElem 输入计算实际单元的函数值时需要的导数信息总个数
 * @param MultiIndexInRefElem 输出需要的RefElem中的导数信息
 * @param IndMultiIndexInRefElem 输出RefElem中的导数信息中 0，1，2阶导数的个数和总个数
 */
void ComputeMultiindexInRefElem(MULTIINDEX *MultiIndexInElem, INT NMultiIndexInElem, 
                                MULTIINDEX **MultiIndexInRefElem, INT *IndMultiIndexInRefElem, INT worlddim);
/**
 * @brief 根据RefElem中所有的函数导数在积分点的值来得到Elem中需要的值
 * @param Base 基函数
 * @param MultiIndex 要计算的微分
 * @param Elem 单元
 * @param BaseValuesInPointsOfRefElem 输入RefElem中所有的函数导数在积分点的值
 * @param IndMultiIndexInRefElem 输入RefElem中的导数信息中 1，2，3阶导数的个数和总个数
 * @param Values 输出Elem中需要的值
 */
void GetBaseValuesWithInfo(BASE *Base, MULTIINDEX MultiIndex, ELEMENT *Elem, DOUBLE *BaseValuesInPointsOfRefElem, 
                            INT *IndMultiIndexInRefElem, DOUBLE *Values);

// 增加边条件(Neumann or Robin)
void AddBoundFun(DISCRETEFORM *DiscreteForm, FUNCTIONVEC *BoundaryFunction, BOUNDARYTYPE bdtype, QUADRATURE *PartialQuadrature);
#endif